Sum-of-Squares matrices

Adapted from: Examples 3.77 of [BPT12]

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R. Semidefinite Optimization and Convex Algebraic Geometry. Society for Industrial and Applied Mathematics, 2012.

Introduction

Consider the symmetric polynomial matrix $P(x) = \begin{bmatrix} x^2 - 2x + 2 & x\\ x & x^2 \end{bmatrix}.$ We could like to know whether $P(x)$ is positive semidefinite for all $x \in \mathbb{R}$.

using DynamicPolynomials
@polyvar x
P = [x^2 - 2x + 2 x
            x     x^2]
2×2 Matrix{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}}:
 2 - 2x + x²  x
 x            x²

A sufficient condition for a symmetric polynomial matrix $P(x)$ to be positive semidefinite for all $x$ is to the existence of a matrix $M(x)$ such that $P(x) = M^\top(x) M(x)$. If such matrix $M$ exists, we say that the matrix is an \emph{sos matrix} (see [Definition 3.76, BPT13]). While determining whether $P(x)$ is positive semidefinite for all $x$, is NP-hard (checking nonnegativity of a polynomial is reduced to this problem for $1 \times 1$ matrices), checking whether $P(x)$ is an sos matrix is an sos program.

using SumOfSquares

We first need to pick an SDP solver, see here for a list of the available choices.

import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

model = SOSModel(solver)
mat_cref = @constraint(model, P in PSDCone())
optimize!(model)
termination_status(model)
OPTIMAL::TerminationStatusCode = 1

While the reformulation of sos matrix to sos polynomial is rather simple, as explained in the "Sum-of-Squares reformulation" section below, there is a technical subtelty about the Newton polytope that if not handled correctly may result in an SDP of large size with bad numerical behavior. For this reason, it is recommended to model sos matrix constraints as such as will be shown in this notebook and not do the formulation manually unless there is a specific reason to do so.

As we can verify as follows, only 3 monomials are used using the sos matrix constraint.

certificate_monomials(mat_cref)
3-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
 ##15454
 x##15456
 x##15454

Sum-of-Squares reformulation

One way to obtain the reduction to an sos program is to create an intermediate vector of variable $y$ and check whether the polynomial $p(x, y) = y^\top P(x) y$ is sos. However, special care is required when approximating the Newton polytope of $p(x, y)$. Indeed, for instance if the entries of $P(x)$ are quadratic forms then the Newton polytope of $p(x, y)$ is the cartesian product between the Newton polytope of $y^\top y$ and the Newton polytope of $x^\top x$. In other words, $p(x, y)$ belongs to a family of quartic forms called biquadratic forms. This fact is important when generating the semidefinite program so that only bilinear monomials are used. So if the cheap outer approximation is used (instead of the exact polyhedral computation) for the newton polytope then it is important to use a multipartite computation approximation of the newton polytope. The multipartie exact approach may perform worse compared to the unipartite exact in certain cases though. Consider for instance the polynomial matrix $\mathrm{Diag}(x_1^1, x_2^2)$ for which $p(x, y) = x_1^2y_1^2 + x_2^2y_2^2$. For this polynomial, only the monomials $x_1y_1$ and $x_2y_2$ are needed in the SDP reformulation while the multipartite approach, as it will compute the Newton polytope as a cartesian product, will not see the dependence between $x$ and $y$ in the presence of monomials and will also select the monomials $x_1y_2$ and $x_2y_1$.

@polyvar y[1:2]
p = vec(y)' * P * vec(y)

\[ 2y_{1}^{2} + 2xy_{1}y_{2} - 2xy_{1}^{2} + x^{2}y_{2}^{2} + x^{2}y_{1}^{2} \]

We can see above that p is biquadratic polynomial in the variables x and y. Computing the Newton polytope with the cheap outer approximation without exploiting this multipartite structure gives the following 6 monomials.

X = monomials(p)
unipartite = Certificate.NewtonDegreeBounds(tuple())
Certificate.monomials_half_newton_polytope(X, unipartite)
6-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
 y₂
 y₁
 x
 y₁y₂
 xy₂
 xy₁

Exploiting the multipartite structure gives 4 monomials.

multipartite = Certificate.NewtonDegreeBounds(([x], y))
Certificate.monomials_half_newton_polytope(X, multipartite)
4-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
 y₂
 y₁
 xy₂
 xy₁

In the example above, there were only 3 monomials, where does the difference come from ? Using the monomial basis, the only product of two monomials that is equal to y[2]^2 is y[2] * y[2]. As y[2]^2 is not a monomial of p, we can conclude that the diagonal entry with row and column corresponding to y[2] will be zero hence the whole column and row will be zero as well. Therefore, we can remove this monomial.

Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(multipartite))
3-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
 y₁
 xy₂
 xy₁

The same reasoning can be used for monomials y[1]y[2] and x therefore whether we exploit the multipartite structure or not, we get only 3 monomials thanks to this post filter.

Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(unipartite))
3-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
 y₁
 xy₂
 xy₁

This page was generated using Literate.jl.